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Summary 

A rapid computational method has been devel- 
oped to determine the equilibrium chemical com- 
position of air to 15 000 K. Eleven chemically re- 
acting species— O 2 , N 2 , O, NO, N, NO + , e~, N + , 
0 + , Ar, and Ar+ — are included. The method in- 
volves the algebraic combination of the equilibrium, 
elemental mass-balance, and charge-neutrality equa- 
tions and solves the resulting equation by an itera- 
tive method. Computational speeds for the determi- 
nation of the equilibrium chemical composition are 
significantly faster than those of the often-used free- 
energy-minimization procedure. Data are also in- 
cluded from which the thermodynamic properties of 
air can be computed. Additionally, a listing of the 
computer program and a set of sample results are 
included. 

Introduction 

The application of computational fluid dynamics 
(CFD) techniques to hypersonic atmospheric flight 
vehicles requires methods that rapidly compute the 
composition of air over a wide range of tempera- 
tures and densities. Most of the available meth- 
ods are very general in nature and employ the free- 
energy-minimization procedure. The generality of 
any method may be desirable if only a few calcu- 
lations are to be carried out. But in a CFD code 
for hypersonic-flow-field computations for which the 
chemistry code is a subroutine executed a very large 
number of times, a large fraction of the computa- 
tional time could be spent in executing the chem- 
istry subroutine alone unless this subroutine is very 
efficient. For such an application it is evident that it 
would be beneficial to forego a certain amount of gen- 
erality in favor of computational speed. The overall 
gain in the computational time, of course, depends 
on the particular application. 

A method of combining the chemical equilibrium 
equations and the mass constraint equations alge- 
braically and solving the resulting equation numeri- 
cally was suggested by Erickson, Kemper, and Alli- 
son (ref. 1). Recently, Erickson and Prabhu (ref. 2) 
applied this method to hydrocarbon combustion 
products. Smith, Erickson, and Eastwood (ref. 3) 
also used a similar method for high-temperature 
equilibrium air. The present method, like that used 
by Smith et al., divides the temperature-density do- 
main into four regions, each with a limited set of 
major species and a different computational scheme. 
Although the present work follows the same general 
approach used by Smith et ah, it differs in that a dif- 
ferent set and a smaller number of major species are 
selected to form the primary iteration equations, and 


all the minor species are included in a series of sec- 
ondary computations. The present approach leads to 
a more rapid computation and ensures continuity of 
results between adjacent regions. Eleven species are 
included, and air is assumed to be a reacting mixture 
of ideal gases. 

The following sections describe the four models, 
their governing equations, and the solution proce- 
dures. A computer program to determine the equilib- 
rium composition of air over a wide range of temper- 
ature and density conditions has been written with 
the present scheme. A listing of this code and a set 
of results are included. Computational speed for the 
present program is significantly faster than that of 
conventional free-energy-minimization routines. 

Symbols 

a^j coefficients of the approximating 

polynomials for thermodynamic 
functions 

b n coefficients of the polynomials 

c p,i specific heat at constant pressure of 

zth species, J/mol-K 

Ef molar internal energy of ith species, 

J/mol 

e specific internal energy of the mixture, 

J/kg 

G® molar Gibbs energy of ith species, 

J/mol 

Hf molar enthalpy of ith species, J/mol 

h specific enthalpy of the mixture, J/kg 

Kj equilibrium constants in terms of mole 

numbers (j = 1, 2, ..., 7) 

K p j equilibrium constants in terms of 

partial pressure (j = 1, 2, ..., 7) 

N(e) number of iterations for an accuracy e 

p pressure of the mixture, Pa 

p° standard-state pressure, 101 325 Pa 

R universal gas constant, 8.314 J/mol-K 

Sf molar entropy of ith species, J/mol-K 

s specific entropy of the mixture, 

J/kg-K 

T temperature of the mixture, K 

T re f reference temperature for enthalpy, 

OK 

t computational time, -1- tiN(e), ps 


t{ time per iteration, ps 

to constant time outside the iteration 

loop, fJ , s 

AG° Gibbs energy for a reaction 

A H® i enthalpy of formation of ith species at 
T re t, J/mol 

AN increase in mole number in a reaction 

e accuracy in mole fraction 

p density of the mixture, kg/m 3 

o total mole number of the mixture, 

mol/kg 

o t mole number of *th species 

(* = 1,2,..., 11), mol/kg 

(7 o , <7 n • 17 A equivalent mole numbers of oxygen, 
nitrogen, and argon, mol/kg 

Analysis 

In the present analysis, low-temperature air is as- 
sumed to consist of O 2 , N 2 , and Ar in molar pro- 
portions of 20.96 : 78.11 : 0.93. Eleven chemically 
reacting species are assumed to be present in equilib- 
rium in air — O 2 , N 2 , O, NO, N, NO + , e~, N + , 0 + , 
Ar, and Ar + — and are designated by the respective 
indices i = 1,2, ..., 11. 

Seven independent chemical reactions involving 
the 11 reacting chemical species are considered. 
These are as follows: 

( 1 ) 0 2 = 20 

(2) N 2 + 0 2 = 2NO 

(3) N 2 = 2N 

(4) NO + + e" = NO 

(5) N+ + e" = N 

(6) 0 + + e" = O 

(7) Ar + + e“ = Ar 

The corresponding 7 nonlinear equilibrium equa- 
tions involving the unknown mole numbers of the 
11 species and the 7 equilibrium constants are the 
following: 


Ki = <73 M 

(1) 

K2 = o\ja \02 

(2) 

b 

II 

£ 

(3) 

K 4 = <74/06 <77 

(4) 


#5 = <75/08(77 

(5) 

Kq = O3/09O7 

(6) 

Ki = <710/(71107 

(7) 


The equilibrium constants Kj in terms of mole num- 
bers are related to the equilibrium constants K p j in 
terms of partial pressures by the following relation: 

Kj = K pJ (p°/pRT) AN 

where AiV is the increase in the number of moles 
when the reaction proceeds from reactants to prod- 
ucts as written for reactions (1) to (7). Note that 

AN = 1 for reactions (1) and (3) 

= —1 for reactions (4), (5), (6), and (7) 

= 0 for reaction (2) 

Also, Kj has the units of moles per kilogram raised 
to the AN power, and K p j has units of atmospheres 
raised to AN power. The equilibrium constants 
K\ to K-j are determined through the use of Gibbs 
energies of formation for each of the chemical species 
participating in the chemical reactions. 

The mass-balance equations for oxygen, nitrogen, 
and argon provide three linear equations involving 
the unknown values of <x, (i = 1 , 2 ,..., 11) and the 


known values of on, 00 , and <r\. 

2<Ji + 03 + 04 + 06 + 09 = Oq (8) 

2(72 + <74 + <75 + <76 + ^8 = (9) 

<710 + <7n = <7\ (10) 

In addition, the charge-neutrality condition of the 
mixture provides the following linear equation: 

<76 + os + < 79 -I- <7ii = <77 (11) 


In principle, it is possible to combine equa- 
tions (1) to (11) into a single equation in one un- 
known and solve by some iteration procedure. The 
basic approach employed in the present analysis, 
however, is to identify 4 regions in the temperature- 
density domain such that in each of these re- 
gions, only some of the 11 species are found in 
abundance. These species are referred to as the 
major species and the rest as the minor species 
in that region. The nonlinear equations govern- 
ing the major species and the appropriate lin- 
ear equations are combined into a single equa- 
tion, with the assumption that the minor species 
are adjustable constants. Such governing equations 
are solved by an iterative procedure updating the 
minor species at each iteration step. (This approach 
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is essentially a combination of Newton and fixed- 
point iteration schemes.) Although some of the 
species are treated as the major species and the oth- 
ers as the minor species, these models for the 4 re- 
gions compute the correct values of all the 11 react- 
ing chemical species considered. These models also 
have some overlapping regions of applicability, and 
the choice of one model over the other for a given 
condition is entirely based on the criterion of mini- 
mum computational time. 

The following sections discuss these four models. 
In each case the governing equations and the solution 
procedures are described. 


5. Compute the trace species: 

&7 = + (^3/Ah) + {V10/K7) 

<76 = <7 4 /(tf 4 <7 7 ) 

<78 = <75/(K5<77) 

<79 = <T3/(Kea 7 ) 

<711 - <y\o/{K 7 (j 7 ) 

As one might expect, this iteration scheme con- 
verges in a couple of iterations at low temperatures, 
where the minor species are present in extremely 
small quantities. 


Model I 

At very low temperatures, the only major species 
are O2, N2, and Ar. Among the remaining species, 
those with a neutral charge (i.e., NO, N, and O) are 
small and are treated as minor species; the ionized 
species (i.e., NO+, e”, N + , 0 + , and Ar + ) are in very 
small quantities and hence are denoted trace species. 
Model I is a very simple model that considers this 
case. The mole numbers of the minor species and 
of the trace species are assumed to be zero to start 
with, and the major species, namely O2, N2, and 
Ar, are computed from the mass-balance equations. 
The equilibrium equations then determine the minor 
species, which are used to update the major species. 
The trace species are computed after this process has 
converged. The steps involved in this simple iteration 
scheme are given in the following: 

1. Assume cr[ = ctq/ 2, a f 2 = <7^/2, and <tio = <ta* 

2. Compute the minor species: 



°t> = \fK3v2 

04 = K e 0 *, 0 Z 


Model II 

This model is applicable for temperatures which 
are high enough for dissociation but not high enough 
for significant ionization. Under this condition, air is 
assumed to be a mixture of O2, N2, O, NO, N, and 
Ar, with the ionized species present only in small 
quantities. The derivation of the governing equation 
and the method employed in solving it are as follows. 
The equations governing the six major species are 
equations (1) to (3) and (8) to (10). With the mole 
numbers of the ionized species ctq to <79 and <7n set 
equal to pseudo constants (equal to zero at the start) , 
these equations are combined into a single equation 
involving one unknown, namely <73. The steps are as 
follows. Equations (1), (2), and (3) yield 


(74/(73(75 = \ JK 2 / KiKs = K e 

Therefore, £74 = K e oz 0 §- Equations (1) and (3) give 

ar = al/Kx 
£ 7 2 = ol/Kz 

Substituting the above into equation (8) and solving 
for (75 in terms of (73 gives 


3. Update the major species: 


<75 = o'o ~ - (73 /K e az 


°l = ( a O - (73 - 04)/% Finally, upon substituting for (74, 02, and (75 in 

(72 = (ctn -04- (75) /2 equation (9) and simplifying, we obtain the following 

quartic equation with <73 as the unknown: 


4. Check and go to step 5 if the results have con- 
verged. If they have not converged, set cr\ = <J\ 
and £7^ = (72 and go to step 2. 
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Kol = 0 
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where 

b 4 = (8 - 2*2)/#! 

63 = 8 - K 2 - 2K 3 K e 

&2 = — 1 8 <tq + 2 K\ + (^o — — K\K%K e 

bx^i-l + KsK^KxJo 
and 

to = 2*1 


7. Compute the minor species: 

<71 = \/ (a 4 /^4) + (cr 5 /i^ 5 ) + (< 73 /^ 6 ) + (ctio/K 7 ) 
< 7 q = a^/K^d'j 
<78 = <7b/Kb<7l 

<79 = <73/ Kq (77 

<711 = <7io! Kiaj 
<710 = <7A - <711 


Note that ct'q = <jq - a§ — ag and a^ = a^ — <tq — 
erg. The equation for 0-3 is solved by the Newton 
iteration scheme. The initial value for the unknown 
0-3 is determined through solution of the following 
equations: 

2 <ti + cr 3 + <74 = <JQ 

ai = vi/Ki 

04/ a 3 = \J 2Ki/a^K2 


8. Check and stop if the results are converged. If 
they are not converged, go to step 3. 

As one might expect, the range of applicability of 
model II overlaps that of model I at the lower temper- 
atures (at about 800 K). Similarly, the applicability 
of this model at higher temperatures overlaps that of 
model III, which is discussed in the next section. 


Model III 


where 04 is an approximate value of 04 . These 
equations imply that all the oxygen is in the species 
O2, O, and NO, which is a good approximation for 
the initial guess in this model. 

The following steps are taken in the solution 
procedure: 

1. Preset ag = ag = ag = 0 and am — 0 a* 

2. Compute a3. 

3. Determine cr' 0 — <jq - ~ ag 

and 0^ — a N “ ~ a S- 

4. Determine b n (n = 0, 1 , . . . , 4). 

5. Compute Aa3 = — F/F 9 , where 

4 

F=J2 b n ° 3 

n = 0 

F’=J 2 n hn^r 1 

n = 1 

6. Update a3 and compute the major species: 


This model is applicable under certain tempera- 
ture and density conditions when oxygen is almost 
completely dissociated and not much ionization has 
taken place. The major species in this case are O, 
N2, N, and Ar, and the remaining seven species 
are treated as minor species. The equations gov- 
erning the major species (eqs. (3) and (8) to (10)) 
are combined into a single equation in one unknown, 
namely a$. This turns out to be a quadratic equation 
for which the coefficients involve the unknown mole 
numbers of the minor species. As in the previous 
case, these are treated as pseudo constants while the 
quadratic equation is solved. After the quadratic is 
solved for as, the mole numbers of the minor species 
are determined and the coefficients of the quadratic 
equation are updated. This process is continued un- 
til a converged solution is obtained. The following 
steps describe this solution procedure: 

1. Assume a f 3 = oq and am = 0'a and set 
04 — ag = as — 0 initially. 

2. Compute a^ = a^ — 04 — ag — ag. 

3. Determine the coefficients b n (n = 0, 1, and 2): 


<73 = 0-3 + Aa3 

<?5 = <70 - ( 2 <r%/ K l) - < 73 ] /Ke<?3 

<74 = #e <73*75 

<72 = (<7n - <74 - <7 5 )/2 

<71 = (<7Q - <74 - <73)/2 


b 0 = -#3^/2 
h = K 3 /2 
b 2 = 1 

4. Solve the quadratic equation 62*75 +64(75 + 60 = 0 
for <75. 
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5. Determine the unknown mole numbers: 

<72 = (<TN - <T5)/2 
<71 = <t'z 2 /Ki 

a 4 = K e (J3(J5 

<77 = y/Wz/Ks) + {(74/ K 4) + (o Fj /K 5 ) + (< 710 JKt) 
(76 — 04I K. 4 G r ] 

<78 = O^jK^cy 

^9 — V3/K6V7 

0\l= ^10 / K'jcr'j 
&10 = &A ~ V\\ 

Finally, 

<73 — “ 2<7i — <74 — <70 — <79 . 

6. Check and stop if results are converged. If they 
are not converged, set 03 = 03 and go to step 2. 

Model IV 

This model considers the very-high-temperature 
case, which is dominated by atomic and ionic species. 
The major species in this model are O, N, e“, N+, 
O+, and Ar, and the others are treated as minor 
species. When equations (5) to (11) governing the 
major species are combined, they reduce to a cubic 
equation involving one unknown, namely <77. The 
minor species mole numbers that appear in the coef- 
ficients of this cubic equation are treated as pseudo 
constants. Starting with a value of zero for the mole 
numbers of the minor species, we solve the cubic 
equation. Next the mole numbers of the remaining 
major and of all the minor species are computed, and 
the coefficients of the cubic equation are updated. 
This iteration process is continued until the results 
are converged. The following steps describe the iter- 
ation procedure: 

1. Preset cri = o<i = 04 = <tq = an = 0. 

2. Determine a f 0 = — 2 a 1 — a 4 — a$ and 

= a N “ 2a 2 G4 ^6- 

3. Determine b n (n = 0, 1, . . . , 3): 

fy) = - v*o~ ^ll “ a 6 

b\ = 1 — — (an + a$)(Kz + K§) 

h = K§ + Kq - (an + g§)KsK§ 

63 - K b Ks 


4. Compute A<77 = — FjF\ where 

F=J2 h *° 7 

n— 0 

F' = nbntf- 1 

71 = 1 

5. Update ay and compute the major species: 

ay = ay + A ay 

<73 = <j' 0 a 7 K 6 /(l + Kq<t 7 ) 

<j 5 = a^a 7 K 5 /{l + K 5 a 7 ) 
o-ll = <7 A /(1 + K 7 a 7 ) 

<710 = o-a - 
<79 = oz/Kq<jq 
<78 = <75/ ^Q><77 

6. Compute the minor species: 

o\ = dl/Ki 
<72 = <75/^1 

<74 — K e asa$ 

^6 = a±/K^ ay 

7. Check and stop if results are converged. If they 
are not converged, go to step 2. 

Table I gives a summary of the major and mi- 
nor species in the four temperature ranges consid- 
ered here. As pointed out earlier, although these four 
models are developed assuming a certain set of major 
species in each, they have overlapping areas of appli- 
cability, and the choice of one model over the other 
for a given temperature and density is based strictly 
on minimum computational time. Figure 1 shows the 
regions in the temperature-density domain, where 
the four models are applicable and give results in 
minimum time. 

Thermodynamic Data 

Computation of thermodynamic properties for 
each of the 11 species as well as for the mixture 
has been included in the program. The data for 
the determination of the thermodynamic properties 
were obtained from several sources. The data for 
the molecular species were obtained directly from 
the tables of reference 4. The heat capacities for 
the monatomic species were computed by Bonnie J. 
McBride of NASA Lewis Research Center using the 
appropriate statistical-thermodynamic relations and 
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the data from references 5, 6 , and 7. The specific 
heat data for all the species were processed, also by 
McBride, using the least-squares method of curve 
fitting employing the approximating polynomials of 
the following forms: 

5 

Cp,i/R = T, a iJ Tj ~ 1 

3 = 1 
5 

Hf/(RT) = £ aijT^/j + a i6 /T 

3=1 

5 

Sf/R = £ aijT^/ij - 1) + a {J In T + 
j= 2 

G°i/{RT) = Hf/(RT) - Sf/R 

The coefficients a z? g and 7 ( i = 1,2, ...,11) are 
the integration constants. The valnes of were 
determined from the enthalpy of formation (taken 
from the JANAF Thermochemical Tables published 
by Dow Chemical Co.) for all the species at the refer- 
ence temperature of 0 K. The values of j were de- 
termined from the values of entropy at T = 298.15 K. 

The coefficients to ( i = 1 , 2 ,..., 11 ) 

are given in table II. The entire temperature range 
from 200 to 15000 K is divided into five smaller 
ranges; these ranges are 200 to 800 K, 800 to 3000 K, 
3000 to 6000 K, 6000 to 10000 K, and 10000 to 
15 000 K. The curves for the approximating least- 
squares polynomials are forced to pass through com- 
mon points at the junctions between the five tem- 
perature ranges, thereby ensuring continuity in the 
values at the junction points. 

The standard-state Gibbs energy for the 11 re- 
acting species is required for the determination of 
the equilibrium constants for the 7 reactions. Equi- 
librium constants K p j for reactions (1) to (7) are 
given by 

K p ,j = exp [—A Gp(RT)] (3 = 1,2,..., 7) 

where AG ( j = v t G® - J2^Gf, i/j, and u[ are the 
i i 

stoichiometric coefficients, and A G® is the standard- 
state Gibbs energy of formation of the zth species. 
For example, for reaction ( 2 ) (N 2 + 02 = 2NO), 

A<?2 = 2(G*°)nc> “ (£°)n 2 “ (^°)o 2 

and 

K p , 2 = exp(-A G° 2 /RT) 


The enthalpy of the mixture in terms of the 
enthalpies of the constituent species and their mole 
numbers is given by 

11 

h = Ufa 

i=l 

Because the reference temperature is 0 K, internal 
energy of each species is simply given by 

E? = H?-RT (t = 1,2 11) 

and the internal energy of the mixture is given by 

1 = 1 

— h — RTo 

The entropy of the mixture is determined from the 
expression 

11 

5 = -Ra\n(p/p°) + Oi [.Sf - R Infer*)] 

1=1 

where p° is the standard-state pressure (101 325 Pa) 
and p is the pressure of the mixture given by 

p — pRTo 

Results and Discussion 

A computer code has been written to implement 
the computational scheme. The thermodynamic data 
in the form of coefficients a^ j (j = 1, 2, . . . , 7) for the 
11 chemical species (i = 1 , 2 , ... , 11 ) are included in 
the code. The program has four subroutines, one for 
each model. This program reads a set of values for 
temperature, density, and convergence criterion and 
computes the chemical equilibrium composition and 
the thermodynamic properties of air. The appendix 
gives a listing of the computer code and a set of 
sample results. 

The total time t for the computation of chemical 
equilibrium composition for any iterative method is 
given by the expression 

t = £0 + 

where to is a constant time required to compute 
the terms outside the iteration loop and includes 
the time for computation of equilibrium constants 
and estimates of initial values to be used for the 


6 



iteration, tj is the time required per iteration, and 
N(e) is the number of iterations required for the 
desired accuracy e. 

A series of cases were run on the CDC 
CYBER 170/750 to determine the three quanti- 
ties to, ti , and N(e) for the present computational 
scheme. Sample results are shown in table III. The 
constant time is essentially the same for all four mod- 
els. However, the time per iteration and the number 
of iterations required for convergence are different for 
the four models. The number of iterations for con- 
vergence depends on the desired accuracy. Over the 
wide range of conditions explored, one or two iter- 
ations were required for a mole fraction accuracy of 
e = 0.001, and in most cases two or three iterations 
were sufficient for an accuracy of e — 0.0001. Re- 
sults for points along the boundaries between various 
models are indicated with an asterisk. 

Computations have been carried out with the 
computer code to obtain equilibrium chemical 


composition for air over a temperature range of 200 
to 15 000 K and for density values of 1 x 10“ 5 , 
1 x 10” 2 , and 1 x 10 2 kg/m 3 . These results are pre- 
sented in figures 2(a), 2(b), and 2(c). 

Concluding Remarks 

The computational scheme presented provides a 
means for rapidly computing the equilibrium chemi- 
cal composition of air to 15 000 K. This scheme was 
developed for the purpose of providing a relatively 
rapid computer code for use with computational fluid 
dynamics programs that deal with high-temperature 
air. A computer code and a thermodynamic data 
base are also included from which the thermody- 
namic properties of air to 15 000 K may be computed. 

NASA Langley Research Center 
Hampton, Virginia 23665-5225 
January 7, 1988 
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Appendix 

Computer Code Listing and Sample Results 


Computer Code 

PROGRAM AIRNEWC INPUT , OUTPUT ,TAPE5= INPUT ,TAPE6=0UTPUT) 
C***************************************************************************'**'*'C 

c 

C THIS PROGRAM COMPUTES THE EQUILIBRIUM CHEMICAL COMPOSITION OF AIR ASSUMING 

C (1) THE DENSITY IS CONSTANT. 

C (2) THE ELEVEN SPECIES: 02, N2, 0, NO, N, N0+, E-, N+, 0+, AR, AND AR+. 

C (3) CURVE FITTED DATA RECEIVED FROM BONNIE J. MCBRIDE. 

C (4) REFERENCE TEMPERATURE IS ZERO KELVIN. 

C ALSO COMPUTES ENTHALPY, INTERNAL ENERGY, AND ENTROPY. 

C 

REAL AACll .7,5) , C( 11 ) , H( 11 ) , EN(11 ) , SC 11 ) , EC 11 ) 

REAL K1 ,K2,K3,K4,K5,K6,K7,KE 

COMMON / A / K1 ,K2,K3,K4,K5,K6,K7,KE 

COMMON / B / SN, SO, SR, EPS 

COMMON / C / S 

COMMON / D / NITER 

EQUIVALENCE ( SI , S( 1 ) ) , ( S2 , S( 2) ) , ( S3 , SC 3 ) ) , ( S4 , S( 4 ) ) , ( S5 , S( 5) ) , 

1 (S6,S(6)) , (S7,S(7)) ,(S8,S(8)) , (S9,S(9)) ,(SA,S(10)) , 

2 (SB, SCll)) 

C 

GBYRT ( I ,K,T)=( AA( I , 1 , K)*( 1 . -ALOG(T) ) -AA( I , 2 , K)*T/2 . -AA( I , 3 , K) 
l+T**2/6.-AA( I , 4 , K)*T**3/12 . 0-AA( I , 5 , K)*T**4/20 . 0+AA( I ,6,K)/T 
1 -AA (I , 7 , K ) ) 

C 

HBYRT (I , K , T ) =AA ( I , 1 , K ) +AA(I,2,K)*T /2 .+AA( I , 3 , K)*T**2/3 . 

1 +AA( I , 4 , K ) *T **3/4 . +AA ( I , 5 , K)*T**4/5 . . 

2 +AA (I , 6 , K ) /T 
C 

SBYR ( I ,K,T)=AA( I ,2,K)*T +AA( I , 3 , K)*T**2/2 ,+AA( I , 4 , K)*T**3/3 . 

1 +AA( I ,5,K)*T**4/4.+AA( I , 1 , K)*ALOG(T)+AA( I ,7,K) 

C 

C THE FOLLOWING CONSTANTS ARE OBTAINED FROM B. J. MCBRIDE, 

C THE BASIC DATA SOURCES ARE: 

C (1) TABULATED DATA FROM: THERMODYNAMIC PROPERTIES OF 

C INDIVIDUAL SUBSTANCES, ACADEMY OF SCIENCES, MOSCOW, USSR, 

C 1978, VOL . 1 , PART 2. (FOR MOLECULAR SPECIES, 

C AND FORMATION VALUES OF ALL SPECIES.) 

C (2) VALUES COMPUTED (BY MCBRIDE) USING STATISTICAL DATA 

C FOR MONATOMIC SPECIES. 

C SOURCE OF BASIC STATISTICAL DATA: 

C (A) NSRDS NBS VOLUME 35, 1971 

C (B) NBS 3, SECTION 5, 1975 

C (C) NBS 3, SECTION 7, 1976 

C 

C NEUTRAL MOLECULAR OXYGEN 02. 

DATA ( (AA( 1 , J , K) , J=1 ,7) , K=1 ,5)/ 

* , 37703733E+01 , - . 2895220 6E-0 2 , . 95332231E-05 , - . 92469913E-08 , 

* . 30191908E-11 , - . 18859751E+Q2 , . 36933495E+01 , 

* . 28969178E+01 , . 23736545E-02 , - . 14917097E-05 , . 46603392E-09 , 

★- . 539451 66E-1 3 , . 82240429E+02 , . 75019384E+01 , 

* . 28421112E+01 , . 13320561E-02 , - . 3391 5853E-06 , . 446521 78E-10 , 

*- . 22914824E-14 , . 58350386E+03 , . 86255038E+01 , 

* . 56821087E+01 , - . 75300588E-03 , . 22980075E-06 , - . 23955920E-10 , 

* . 80048472E-15 , - . 24700372E+04 , - . 98740 054E+01 , 
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12025200E-10 , 


I 

I 


*- . 27258966E+00 , . 20115141E-02 , - . 24547714E-06 , . 

*— . 21389359E-15 , . 76117740E+04 , . 31631 739E+02/ 

C NEUTRAL MOLECULAR NITROGEN N2 

DATA ( (AA(2, J,K) , J=1 ,7) ,K=1 ,5)/ 

* . 34622648E+01 , . 58202350E-03 , - . 30525456E-05 , . 62280057E-08 , 

*- . 33755957E-11 , . 87951739E+00 . . 32192648E+01 , 

* . 27022403E+01 , . 19443933E-02 , - . 89300044E-06 , . 19739206E-09 , 
. 1 696781 3E-1 3 , . 20180222E+03 , . 72040846E+01 , 

* . 391 4350 6E+01 , . 31537108E-03 , - . 56481043E-07 , . 36012299E-11 , 

* . 52359433E-16 , - . 53551958E+03 , . 21671998E-01 , 

* . 12657466E+01 , . 18269791E-02 , - . 37583896E-06 , . 33033716E-10 , 
*- . 93651108E-15 , . 31426960E+04 , . 17943292E+02 , 

* . 27715943E+02 , - . 74173486E-02 , . 82395970E-06 , - . 35285300E-10 , 

* . 49671195E-15 , - . 56942692E+05 , - . 17402873E+03/ 

C NEUTRAL MONATOMIC OXYGEN 0 

I DATA ( (AA<3, J,K) , J=1 ,7) ,K=1 ,5)/ 

* . 32167143E+01 , - . 37822689E-02 , . 84746788E-05 , - . 88658246E-08 , 

* . 35365597E-11 , . 29640462E+05 , . 1852641 5E+01 , 

* . 26045368E+01 , - . 17235464E-03 , . 115741 38E-06 3641 7858E-10 , 

* . 46011536E-14 , . 29728962E+05 , . 45865262E+01 , 

* . 281 01683E+01 , - . 29039910E-03 , . 90833347E-07 , - . 99427817E-11 , 

* . 3770 4361 E-l 5 , . 29536612E+05 , . 32536467E+01 , 

* . 19209266E+01 , . 21776554E-03 , - . 18288405E-07 , . 49050959E-12 , 

* . 28507296E-17 , . 30783423E+05 , . 92748629E+01 , 

* . 19209266E+01 , . 21776554E-03 , - . 18288405E-07 , . 49050959E-12 , 

* . 28507296E-17 , . 30783423E+05 , . 92748629E+01/ 

C NEUTRAL NITRIC OXIDE, NO 

DATA ( (AA(4, J,K) , J=1 .7) ,K=1 ,5)/ 

* . 420 64360 E+01 , - . 45098354E-02 , . 10557385E-04 , - . 85919382E-08 , 

* . 24047097E-11 , . 10888965E+05 , . 23137934E+01 , 

* . 27543778E+01 , . 23093284E-02 , - . 12823359E-05 , . 34043524E-09 , 

*- . 34807544E-13 , . 11134325E+05 , . 90789671E+01 , 

* . 38015418E+01 , . 49857546E-03 , - . 12531319E-06 , . 14093893E-10 , 

. 44820194E-15 , . 10666566E+05 , . 31619390E+01 , 

* . 49133164E+01 , . 61755267E-06 , - . 55222382E-07 , . 11686489E-10 , 

*- . 54642492E-15 , . 88453776E+04 , - . 45786900E+01 , 

* . 20456654E+02 , - . 61498079E-02 , . 86914275E-06 , - . 50877512E-10 , 

* . 10624157E-14 , - . 22955242E+05 , - . 11561956E+03/ 

C NEUTRAL MONATOMIC NITROGEN N 

DATA ( ( AA ( 5 , J , K ) , J=1 ,7) ,K=1,5)/ 

* . 250 00 000 E+01 , 0 . ,0. ,0. 

*0. , . 56626706E+05, . 41807280E+01 , 

* . 25075109E+01 , - . 24797882E-04 , . 23641517E-07 , - . 1 5288108E-10 , 

* . 28913714E-14 , . 56624949E+05 , . 41431864E+01 , 

* . 263760 51 E+01 , - . 87373323E-05 , - . 64772651 E-0 7 , . 23473436E-10 , 

*— . 1 73961 64E-1 4 , . 56452267E+05 , . 32232104E+01 , 

* . 33720620E+01 , - . 88554639E-03 , . 25293270E-06 , - . 23187896E-10 , 

* . 7047142QE-15, . 56270151E+05 , - . 10563976E+01 , 

*- . 10205642E+02 , . 42931 347E-02 49310664E-06 , . 24951067E-10 , 

★- . 46969943E-15 , . 84933779E+05 , . 96403744E+02/ 

C UNI POSITIVE NITRIC OXIDE, NO+ 

DATA ( ( AA ( 6 , J , K ) ,J=1,7) ,K=1,5)/ 

* . 35659323E+01 27420 30 7E-03,-. 52284273E-06 , . 30228201E-08 , 

191 1270 2E-11 , . 11841356E+06 , . 35935691E+01 , 

* . 27152625E+01 , . 18998146E-02 , - . 85346216E-06 , . 18390288E-09 . 


9 


+-.15323469E-13, , 118621 29E+06 , . 793691 46E+0 1 , 

* . 39310989E+G1 , . 30660785E-03 , - . 567G5942E-07 , . 45573508E-11 , 

+-.77762022E-16, .11786296E+06, . 70198182E+00 , 

* . 1 06561 63E+01 , . 17917961E-02 , - . 33077778E-G6 , . 25052719E-10 , 

- 54540 127E -15 , . 12214242E+06 , . 20328302E+G2 , 

* . 31376951 E+02 , - . 78170663E-02 , . 754371 24E-06 , - . 247341 42E-10 , 

* . 1 5952361 E-15, . 48124504E+05 , - . 20218492E+03/ 


C ELECTRON, E- 

DATA ( (AAC7, J,K) , J=1 ,7) ,K=1 ,5)/ 

* .25GG0000E+01 ,0 . ,0. ,0. 

*G . , - . 42498957E-Q1 1173391 7E+Q2, 

* . 250 00000 E+0 1,0. ,0. ,0. 

★0 . , - . 42493957E-01 , - . 1173391 7E+02 , 

* .25000000E+Q1 ,0 . ,0. ,0. 

*0 . ,-.42498957E-01 , - . 11733917E+G2 , 

* . 25G00G00E+01 , 0 . ,0. ,0. 

*0 . ,-.42498957E-01 , 11 73391 7E+02 , 

* , 2500 00 00 E+0 1 , 0 . ,0. ,0. 

*0 . , - ■ 42498957E-01 , - . 11733917E+02/ 

C UNI POSITIVE MONATOMIC NITROGEN, N+ 


DATA ( (AA(8, J,K) , J=1 ,7) ,K=1 ,5)/ 

* . 30816977E+01 , - . 37737512E-02 , . 96771 466E-0 5 , - . 11099460E-07 , 

* .47248633E-11 , . 2253280 7E+0 6 , . 2431 5936E+0 1 , 

* . 253470 05E+01 , - . 55067384E-04 , . 364400 13E-07 1231 31 99E-10 , 

* , 19224650 E-l 4 , . 22539524E+06 , . 47895586E+0 1 , 

* . 280 05551 E+0 1 , - . 29566365E-03 , . 96649663E-07 11 053590 E-10 , 

* . 441 451 51 E-l 5 , . 2251 8494E+0 6 , . 31 30 5296E+0 1 , 

* . 20047667E+01 , . 191 89730 E-03 , - . 154991 36E-07 , . 41165940E-12 , 

* , 26287640 E-l 7 , . 22622600E+06 , . 84634959E+0 1 , 

* . 20598363E+01 , . 1 8875424E-03 , - . 1S075532E-07 , . 80284344E-12 , 

. 1 3089587E-16 , . 22602766E+06 , . S0254371E+01/ 

C UNI POSITIVE MONATOMIC OXYGEN, 0+ 

DATA ( ( AA ( 9 , J , K ) , J=1 ,7) ,K=1 ,5)/ 

* . 25GG0000E+01 , 0 . ,0. ,0. 

*0 . , .18771 321 E+0 6 , . 4380 1 998E+0 1 , 

* . 25013354E+01 , - . 40431297E-05 , . 43271 450E-08 19527702E-11 , 

+ . 31632791E-15, . 18771288E+06 , . 43734240 E+0 1 , 

* . 22258863E+01 , . 296501 43E-03 11 5611 44E-06 , . 18461 192E-10 , 

+-.89249152E-15, . 1 8791 1 59E+0 6 , . 60 5761 33E+0 1 , 

* . 41210 991E+01 ,-.867861 75E-0 3, . 15337469E-06 , - . 92260135E-11 , 

* .17342049E-15, . 18543702E+06 , - . 66388864E+01 , 

* . 17349474E+01 , - . 80710732E-04 , . 5980 631 7E-0 7 , - . 45998476E-11 , 

* . 10295182E-15 , . 191 07440E+06 , . 1 0791908E+02/ 

C NEUTRAL ARGON, A. 

DATA ( CAAC10 , J,K) , J=1 ,7) , K=1 ,5)/ 

* .25000000E+01 ,0 . ,0. ,0. 

*0 . , - . 42498957E-Q 1 , . 43664979E+0 1 , 

* . 25000000E+01 ,0 . ,0. ,0. 

★0 . , - . 42498957 E-0 1 , . 43664979E+0 1 , 

* . 25000000E+01 ,0 , ,0. ,0. , 

+0 . , - . 42498957E-0 1 , . 43664979E+0 1 , 

* . 26303545E+01 77393441 E-0 4 , . 1 721 9729E-0 7 , - . 170 28231 E-l 1 , 

* . 63198847E-16 , - . 175480 10 E+0 3 , . 34890 096E+01 , 

+-.91423987E+01 , . 385081 69E-02 46S0044GE-06 , . 2450 691 8E-10 , 

*- . 45646899E-15 , . 2775191 5E+05 , . 89461 765E+02/ 
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C UNI POSITIVE MONATOMIC ARGON, A. 

DATA ( (AACll , J,K) , J=l,7) , K=1 , 5)/ 

* . 26818797E+01 , - . 20 720 300E-02 , . 74786Q02E-05 , - . 8767201 3E-03 , 

* . 34484189E-11 , . 18287G91E+06 , . 50764556E+01 , 

* . 24220 6 46 E+ 01 , . 86346979E-03 , - . 78459395E-06 , . 26673998E-09 , 

*- . 31850590E-13 , . 1S285266E+06 , . 59236280E+01 , 

* .25588416E+01 , . 83736230E-04 50091 586E-07 , , 882981 98E-11 , 

*- . 52371 930 E~ 15 , . 1830138GE+06 , . 55313G33E+01 , 

* .25G14164E+Q1 , . 2944655SE-04 , - . 80 471729E-08 , . 77535785E-12 , 

*- . 25780935E-16 , . 18325162E+Q6 , . 60544022E+Q1 , 

* . 25706058E+01 , - . 20 6210 51 E-Q 4 , . 297860 13E-08 , - . 2110 1899E-1 2 , 

* , 57476758E-17 , . 1 8322322E+06 , . 56165023E+01/ 

RR=S. 31441 

P0=1 .01325E+05 
C 

C ASSUME AIR TO HAVE 02 : N2 : A = 20 . 96 : 78.11 s 0.93 

C 

SO = 14.4802 

SN = 53.9620 

SR = 0.3212 

C******************************************************************************C 

c 

C INPUT DATA: 

C 

C TKEL = TEMPERATURE (KELVIN) 

C RKGM3 = DENSITY (KG/M**3) 

C EPS = CONVERGENCE CRITERION 

C 

C******************************************************************************C 
60 NR I TE( 6 , 66) 

66 FORMAT (/,2X,* INPUT:- TKEL, RKGM3 , EPS*, 

1 / , 2X , * ENTER < RETURN) TO STOP*) 

READ (5,*) TKEL, RKGM3, EPS 
I F( EOF ( 5) ) 62,5,62 
5 CONTINUE 
C 

C CHOOSE APPROPRIATE CONSTANTS DEPENDING UPON THE TEMPERATURE. 

C 

K=1 

I F(TKEL .GT. 800.) K=2 
I F (TKEL .GT. 3000.) K=3 
I F(TKEL .GT. 6000.) K=4 
I F (TKEL .GT. 10000.) K=5 
C 

C ..... COMPUTE DG/RT, H/RT, AND S/R FOR ALL THE ELEVEN SPECIES. 

C 

DO 8 1=1,11 

C( I ) =DGBYRT ( I , K ,TKEL) 

H( I ) =( HBYRT ( I , K ,TKEL)*TKEL)*RR 
EN( I ) =SBYR (I , K , TKEL) *RR 
C NRITE( 6,50) I , K ,TKEL , C( I ) , H( I ) , EN( I ) 

50 FORMAT (212, 7E12 .5) 

8 CONTINUE 

C ******************************************************************************C 

c. 


I 

i 
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noon 


C COMPUTE THE EQUILIBRIUM CONSTANTS FOR THE GIUEN TEMPERATURE 

C AND DENSITY CONDITIONS. 

C.....THE REACTIONS ARE: 

C j 

C 02 = 2 0 (1) ! 

C 

C N2 + 02 = 2 NO (2) 

C 

C N2 = 2 N (3) 

C j 

C ( NO*f ) + E- = NO (4) 

C 

C (N+) + E- = N (5) i 

C < 

C (0+) + E- = 0 (6) I 

C 

C (AR+) + E- = AR (7) 1 

c ! 


Ck*kkkkkkkkkkkkkkkkTkkkkkkkkkkkkkkk*kkkkkkkkk~kk'k’k m k-kkkkkkk’k'kicki<:kk'kick m k m k~kkk:k:kjck-kk'kkC 
Kl-C( 3)+C( 3) -C( 1 ) 

K2=C(4)+C(4)-C(1)-C(2) 

K3=C( 5)+C( 5) -C( 2) 

K4=C(4)-C(6)-C(7) 

K5=C( 5) -C( 8) -C( 7) 

K6=C( 3) — C( 9) ~C( 7) 

K7=C(10)-C(11)-C(7) 

K1=EXXP( -K1 ) * K2=EXXP( -K2) * K3=EXXP(-K3) 

K4=EXXP ( -K4 ) * K5=EXXP( -K5) * K6=EXXP(-K6) 

K7=EXXP( -K7) 

P OBYRT=PO/( RR*TKEL*RKGM3) 

K1=K1*P0BYRT 
K3=K3*PGBYRT 
K4=K4/P0BYRT 
K5=K5/P0 BYRT 
K6=K6/P0BYRT 
K7=K7/P0BYRT 
KE « SQRT ( K2/C K1*K3) ) 

CALL APPROPRIATE SUBROUTINE DEPENDING UPON THE TEMPERATURE 

TO COMPUTE THE COMPOSITION. 

IFCTKEL .LE. 1600.0) GO TO 151 
RLOG = ALOGIO ( RKGM3) 

TEMPI = 1000 .0*C3. 5+10 .0*C C RLOG+6 , 0 )/7 . 0)**2. 5) 

IFCTKEL .LE. TEMPI) GO TO 154 

TEMP2=1000 .*( 5. 4+3. 2*C RLOG+6. )/6. +5. 55*CABS( RLOG+6. )/6. )**3. 35) 

IFCTKEL .LE. TEMP2) GO TO 152 

153 MODEL = 3 
CALL CUBIC 
GO TO 156 

152 MODEL = 2 

CALL QUDRATC 
GO TO 156 

154 MODEL = 4 
CALL QUART I C 
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GO TO 156 
151 MODEL = 1 

CALL LOWTEMP 
C 

C TOTAL MOL NUMBER. 

C 

156 SUM - S1+S4+S2+S3+S5+S6+S7+S8+S9+SA+SB 
C 

C COMPUTE PRESSURE 

C- 

PBAR = RKGM3*8 . 314*SUM*TKEL/100Q0Q . 

PATM = P BAR/1 . 0 1 325 
C 

C COMPUTE RESIDUES.... 

C 

R1 = S3*33/K1/S1 -1 .0 
R2 = S4*S4/K2/S1/S2-1.G 
R3 = S5*S5/K3/S2-1 .0 
R4 = S4/K4/S6/S7-1 .0 
R5 = S5/K5/S8/S7-1 .0 
R6 = S3/K6/S9/S7-1 . 0 
R7 = S0-2,0*S1-S3-S4-S6-S9 
R8 = SN-2 . G*S2-S4-S5“S6-S8 
R9 = S7-S6-S8-S9-SB 
RA = SR - SA -SB 
RB = SA - K7*SB*S7 
I F ( SB .EQ. 0.0) RB = 0.0 

ikirM-k-kQ 

C 

C. ... .COMPUTE THE ENTHALPY, INT. ENERGY, AND ENTROPY OF THE MIXTURE. 

C NOTE: ENTHALPY AND INT. ENERGY ARE IN JOULES/KG. 

C ENTROPY IS IN JOULES/KG . KELU IN 

C 

HO = 0.0 
EO = 0.0 

ENTO = -RR*SUM*ALOG(PATM/SUM) 

DO 96 1=1,11 
HO = H0+S(I)*H(I) 

E(I) * H( I ) -RR*TKEL 
E0=E0+S( I )*E( I ) 

I F ( S ( I ) .EQ. 0.0) GO TO 96 

ENTO = ENTO+SC I )*(EN( I )-RR*ALOGC S( I )) ) 

96 CONTINUE 
C 

C WRITE OUT THE RESULTS. 

C 

WRITEC 6,101) TKEL , RKGM3 , EPS 

WRITE( 6,100) SO , SN , SR 

WRITEC 6, 102) K1 , K.2 , K3 , K4 , K5 , K6 , K7 , KE 

I F (MODEL .EQ. 1) WRITEC6,304) NITER 

I F (MODEL .EQ. 2) WRITE(6,307) NITER 

I F (MODEL .EQ. 3) WRITE(6,306) NITER 

I F (MODEL .EQ. 4) WRITE(6,3G5) NITER 

WRITE (6, 30 3) SI , S2 , S3 , S4 , S5 , S6 , S7 , S3 , S9 , SA , SB 
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WRITE( 6,318) ( H( I ) , I =1,11) 

WRITE( 6 , 316) (E(I), I =1,11) 

WRITE( 6 , 319) ( EN ( I ) , I =1,11) 

WRITEC 6,317) SUM , HO , ENTO , EO , PBAR 

WRI TEC 6,204) R1 , R2 , R3 , R4 , R5 , R6 , R7 , R8 , R9 , RA , RB 

303 FORMAT (/ , 2X , 38HM0LE NUMBERS OF THE SPECIES (MOL/KG):-, 

1/, 9X , 2H02 , 1 IX , 2HN2 , 1 IX , 2H0 , 1 IX , 2HN0 , 1 IX , 2HN ,10X,3HNO+, 

1/ , 9X , 2HE- , 11X , 2HN+ , 11X , 2H0+ , 11X , 2HAR , 11X , 3HAR+ ,// , 

1 2X , 6E13 . 5 ,/ , 2X , 5E13 . 5) 

304 FORMAT (/ , 2X , 43HS0LUTI ON BY THE SIMPLE ITERATION PROCEDURE;, 

1 I 3 , 1 1H ITERATIONS) 

305 F ORMAT ( / , 2X , 33HS0LUT I ON OF THE QUARTIC EQUATION;, 

1 I 3 , 11H ITERATIONS) 

306 FORMAT (/ , 2X , 31HS0LUTI ON OF THE CUBIC EQUATION;, 

1 I 3 , 1 1 H ITERATIONS) 

307 FORMAT (/ , 2X , 35HS0LUTI ON OF THE QUADRATIC EQUATION;, 

1 I 3 , 11H ITERATIONS) 

100 FORMAT (/,2X,*0XYGEN = * , F8 . 4 , 2X ,*NITROGEN = *,F8.4,2X, 

1 JARGON = * , F6 . 4 ,* MOLS/KG.*) 

101 FORMAT (/,2X,*T (KELVIN) = * , F8 . 2 , 2X ,*DENSITY (KG/M**3) = *,E12.6 
1 ,2X,*EPS = *,E10.4) 

204 FORMAT (/ , 2X ,*RESI DUES : , 3(/ , 2X , 6E13 . 5) ) 

102 FORMAT (/,2X,*EQUI LIBRIUM CONSTANTS K(I):-*, 

1 / , 2X , 6E13 . 5 ,/ , 2X , 6E13 . 5) 

317 FORMAT ( 

1/ , 2X , 29HT0TAL MOL NUMBER = ,F12.4,11H MOLS/KG. 

1/ , 2X , 29HENTHALPY OF THE MIXTURE = ,E12.6,11H JOULES/KG., 

1/ , 2X , 29HENTR0PY OF THE MIXTURE = ,E12.6,18H JOULES/KG ./KELVIN , 

1/ , 2X , 29HINT . ENERGY OF THE MIXTURE = ,E12.6,11H JOULES/KG., 

1 / , 2X , 29HPRESSURE = ,E12.6,6H (BAR)) 

318 FORMAT (/ , 2X , 39HENTHALPY OF EACH SPECIES (JOULES/MOL):-,/, 

1 3 ( 2X , 6E1 3 . 5 , / , ) ) 

316 FORMAT (/,2X,46H INTERNAL ENERGY OF EACH SPECIES (JOULES/MOL):-,/, 

1 3 ( 2X , 6E1 3 . 5 , / , ) ) 

319 FORMAT (/ , 2X , 45HENTP.0PY OF EACH SPECIES (JOULES/MOL/KELVIN):-,/, 

1 3 ( 2X , 6E1 3.5,/,)) 

GO TO 60 
62 STOP 
END 
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SUBROUTINE LOWTEMP 

C ************************************************ 


THIS SUBROUTINE COMPUTES THE COMPOSITION OF AIR 
AT LON TEMPERATURES NHERE THERE IS SOME DISSOCIATION 
BUT MERY LITTLE IONIZATION. 

THE MAJOR SPECIES ARE 02, N2, AND AR . 

(MODEL I) 


jcJc&*jcJeTkjcTkjcjcTk*jcjcft&&Tkik*jcjc1cjc*jcjc*jcjc1cjc-k1cic1cjcic*1cicjcjc*jc1cicicjcTk1cjcjcJcjcjcjc*1cjc*icjc1cjcicjc1cicicjcJcjcicjc&-kC 
COMMON / A / K1 ,K2,K3,K4,K5, K6,K7,KE 

COMMON / B / SN, SO, SR, EPS 

COMMON / C/ SI , S2 , S3 , S4 , S5 , S6 , S7 , SB , S3 , SA , SB 

COMMON / D / ITER 

REAL K1 ,K2,K3,K4,K5,K6,K7,KE 


INITIALIZE AND START ITERATING. 


SIP = SO / 2.0 * S2 = SN / 2.0 * SA = SR 

DO 200 ITER =1,20 

COMPUTE THE MINOR SPECIES. 


S3 = SQRT(K1*S1P) * S5 = SQRT ( K3*S2) * S4 = KE*S3*S5 


UPDATE THE MAJOR SPECIES. 


SI = (SO - S4 - S3) / 2.0 * S2 = ( SN - S4 - S5) / 2.0 

CHECK FOR CONVERGENCE. 

DS1 = SI - SIP 

I F ( ABS( DS1 ) .LE. EPS) GO TO 201 
SIP = SI 

200 CONTINUE 

201 CONTINUE 


COMPUTE THE TRACE SPECIES. 


C3 

= S3/K6 

ii 

o 

* 

S4/K4 

* 

C5 = S5/K5 

* 

CG = SA/K7 

S' 7 

= SQRT ( C3 + C4 + C5 

+ C6) 





S6 

= C4/S7 

* SB = 

C5/S7 

* 

S3 = C3/S7 

% 

SB = C6/S7 


RETURN 

END 
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SUBROUTINE QUART I C 

C******************************************************************************C 


THIS SUBROUTINE COMPUTES THE COMPOSITION OF AIR 
AT MODERATE TEMPERATURES WHERE THERE IS DISSOCIATION 
BUT NOT MUCH IONIZATION. 

THE MAJOR SPECIES ARE 02, N2, 0, NO, AND N. 

(MODEL II) 


COMMON / A / K1 , K2 , K3 . K4 , K5 , K6 , K7 , KE 
COMMON / B / SN, SO, SR, EPS 

COMMON / C / SI , S2 , S3 , S4 , S5 , S6 , S7 , S8 , S9 , SA , SB 

COMMON / D / ITER 

REAL K1 ,K2,K3,K4,K5,K6,K7,KE 


INITIALIZE, DETERMINE THE INITIAL GUESS, AND START ITERATION. 

S6 = 0.0 * SB =0.0 * 39 = 0.0 * SA=SR 

B4 = (8.0 - K2 - K2) / K1 

B3 = 8.0 - K2 - 2 . 0*K3*KE 

BB= Kl*( 1.0 + SQRT ( K2*SN / (K1 + Kl))) 

S3 = ( - BB + SQRT ( BB*BB + 8.0*S0*K1)) / 4.0 

C1=K1*K3*KE 

DO 200 ITER = 1,20 

SOI = SO - S6 - S9 

SN1 = SN - S6 - S8 

B2 = - 8 . 0*S01 + Kl + Kl + (SOI - SN1)*K2 - Cl 
B1 = S01*( - Kl - Kl - Kl - Kl + Cl) 

80 = 2.0*S01*S01*K1 

F = BO + S3*(B1 + S3*(B2 + S3*(B3 + S3*B4) ) ) 

FP = B1 + S3*(B2 + B2 + S3*(B3 + B3 + B3 + S3*( B4+B4+B4+B4) ) ) 
D3 = - F / FP 
S3 = S3 + D3 


. UNRAUEL THE 

REST 

OF THE MAJOR SPECIES 

S5 

s- 

(SOI - 

S3 - 

2.G*S3*S3 / Kl) / ( KE*S3) 

S4 

SB 

KE*S5*S3 


S2 

as 

( SN1 - 

S4 - 

S5) / 2.0 

SI 

= 

(SOI - 

S4 - 

S3) / 2.0 


COMPUTE THE MINOR SPECIES. 


C3 

S7 

= 

S3/K6 * C4 = 

SQRT ( C3 + C4 + C5 

S4/K4 
+ C6) 

% 

C5 = S5/K5 

* 

C6 

= SA/K7 

S6 

SA 

s 

C4/S7 * S8 « 

SR - SB 

C5/S7 

% 

S3 = C3/S7 

* 

SB 

= C6/S7 


CHECK FOR CONVERGENCE 

I F( ABS( D3) .LE. EPS) GO TO 201 

200 CONTINUE 

201 CONTINUE 
• RETURN 

END 


I 
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non ono non non oonoooo 


SUBROUTINE CUBIC 

C********tHc***************tHt*tWc**************tHt*******^ 

THIS SUBROUTINE COMPUTES THE COMPOSITION OF AIR 

AT MERY HIGH TEMPERATURES WHERE NOT MANY DIATOMIC PARTICLES ARE FOUND. 
THE MAJOR SPECIES ARE 0, N, E- , N+ , AND 0+ . 

(MODEL IM) 


COMMON / A / K1 , K2 , K3 , K4 , K5 , K6 , K7 , KE 
COMMON / B / SN, SO, SR, EPS 

COMMON / C / SI ,S2,S3,S4,S5,S6,S7,S8,S9,SA,SB 

COMMON / D / ITER 

REAL K1 ,K2,K3,K4,K5,K6,K7,KE 


INITIALIZE, DETERMINE THE INITIAL GUESS, AND START ITERATION. 

SI =0.0 $ S2 = 0.0 * S4 = 0.0 * S6 = 0.0 

SA = SR % SB = 0.0 

S7 = ( SQRTd.O + 4 . 0*SN*K5) -1 . 0 ) / ( K5 + K5) 

S3 = SO / (1.0 + K6*S7) 

SB = SR / (1.0 + K7*S7) 

S7 = S7 + S3 + SB 

B3 = K5*K 6 

C2 = K5 + K6 

DO 200 ITER = 1,20 

SOI = SO - SI - SI - S4 - S6 

SN1 = SN - S2 - S2 - S4 - S6 

S6P = S6 + SB 

B2 = C2 - B3*S6P 

B1 = 1.0 - SNICKS - S01*K5 - S6P*C2 


BO 

= 

- SN1 - SOI 

- S6P 


F 

= 

BO + S7*(B1 

+ S7*(B2 + S7*B3) ) 

FP 

= 

B1 + S7*(B2 

+ B2 + 

S7*3.0*B3) 

D7 

= 

- F / FP 



S7 

as 

S7 + D7 



. . UN RAM EL THE REST 

OF THE 

MAJOR SPECIES 

S3 

= 

S01*S7*K6 / 

(1.0 + 

K6*S7) 

S5 

= 

SN1*S7*K5 / 

(1.0 + 

K5*S7) 

SB 

= 

SR / (1.0 + 

K7*S7) 


SA 

= 

SR - SB 



S9 

= 

S3 / (K6*S7) 



S3 

= 

S5 / (K5*S7) 




....COMPUTE THE MINOR SPECIES. 

51 = S3*S3 / K1 

52 = S5*S5 / K3 
S4 = KE*S3*S5 

S6 = S4 / ( K4*S7) 

CHECK FOR CONMERGENCE 

I F ( ABS ( D7 ) .LE. EPS) GO TO 201 

200 CONTINUE 

201 CONTINUE 
RETURN 
END 
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non non oooo 


SUBROUTINE QUORATE 
C 

C THIS SUBROUTINE COMPUTES THE COMPOSITION OF AIR AT INTERMEDIATE 

C. ... .TEMPERATURES WHERE OXYGEN IS ALMOST FULLY DISSOCIATED, 

C BUT NOT MUCH IONIZATION HAS TAKEN PLACE. 

C THE MAJOR SPECIES ARE N2 , 0, AND N. 

C (MODEL III) 

C 

COMMON / A / K1 ,K2,K3,K4,K5,K6,K7,KE 

COMMON / B / SN, SO, SR, EPS 

COMMON / C/ SI , S2 , S3 , S4 , S5 , S6 , S7 , S8 , S9 , SA , SB 

COMMON / D / ITER 

REAL K1 ,K2,K3,K4,K5,K6,K7,KE 

INITIALIZE, DETERMINE THE INITIAL GUESS, AND START ITERATION. 

54 s= 0.0 * S6 = 0.0 * SB = 0.0 * SA = SR 

B1 = K3 / 2.0 * B2 = K1 / 2.0 

55 = ( - B1 + SQRT(B1*B1 + 4.0*B1*SN )) / 2.0 

S3P = ( - B2 + SQRT(B2*B2 + 4 . 0*B2*SO )) / 2.0 

DO 200 ITER = 1,20 

SN1 = SN - S4 - S6 - S8 
F = - B1*SN1 + S5*(B1 + S5) 

FP = B1 + S5 + S5 
S5 = S5 - F / FP 
S2 = ( SN1 - S5) / 2.0 

UNRAVEL THE REMAINING SPECIES 


SI 

= 

S3P*S3P / 

K1 






S4 

= 

KE*S3P*S5 







C3 

= 

S3P/K6 

* C4 

= S4/K4 

* 

C5 = S5/K5 

% 

06 = SA/K7 

S7 

= 

SQRT ( C3 + 

04 + C5 

+ C6) 





S6 

= 

C4/S7 

* S8 

= C5/S7 

* 

S9 = C3/S7 

* 

SB = C6/S7 

SA 

= 

SR - SB 







S3 

= 

SO - SI - 

SI - S' 4 

- S6 - 

S9 





CHECK FOR CONVERGENCE. 

D3 = S3P - S3 
S3P = S3 

I F( ITER .EQ. 1) GO TO 200 
I F( ABS(D3) .LE. EPS) GO TO 201 

200 CONTINUE 

201 CONTINUE 
RETURN 
END 

FUNCTION EXXP(X) 

EXXP=0 .0 

IFCX .LT. - 670.0) RETURN 
EXXP=EXP(X) 

RETURN 

END 
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Sample Results 

/lqo 


INPUT : - TKEL , RKGM3 , EPS 

ENTER < RETURN) TO STOP 
? 10000 . ,1 .e-6,1 .e-10 

T (KELVIN) = 10000.00 DENSITY ( KG/M 3 ) = .1QQQQQE-Q5 EPS = .100GE 

OXYGEN = 14.4802 NITROGEN = 53.9620 ARGON = .3212 MOLS/KG. 

EQU I L I BR I UN CONSTANTS K( I ) s - 


. 2950 8E4-1 1 

. 1 9149E+G1 

. 23078E+09 

. 34684E-04 

. 13169E-02 

. 19413E-02 

. 53027E-09 




SOLUTION OF THE 

CUBIC EQUATION-, 5 ITERATIONS 


MOLE NUMBERS OF 

THE SPECIES 

(MOL/KG) : - 



02 

N2 

0 

NO 

N 

E- 

N+ 

0+ 

AR 

AR+ 

. 9100 3E-10 

. 73958E-0 7 

. 16387E+01 

. 35900E-08 

. 41 314E+01 

. 62958E+02 

. 49831E+ 02 

. 12841E+02 

. 34983E-01 

. 28622E+0Q 

ENTHALPY OF EACH 

1 SPECIES (JOULES/MOL):- 



. 41091E+06 

.38016E+06 

.46617E+06 

. 48132E+06 

. 71628E+06 

. 2G786E+06 

. 2G934E+G 7 

. 1 7866E+Q7 

. 20790E+06 

. 1 7334E+07 

INTERNAL ENERGY 

OF EACH SPECIES (JOULES/MOL):- 


. 32776E+ 06 

. 29702E+06 

. 38303E+ 06 

. 39817E+06 

. 63313E+06 

. 12472E+06 

. 20103E+07 

. 1 7035E+07 

. 12476E-H06 

. 1 6502E+07 

ENTROPY OF EACH 

SPECIES (JOULES/MOL/KELVIN):- 


. 3362OE+03 

. 31 386E+03 

. 23614E+G3 

. 33631E-H 0 3 

. 231 35E+03 

. 93886E+02 

. 23460E+03 

. 2301 3E+03 

, 22776E+03 

. 24261E+ 03 

TOTAL MOL NUMBER 

- 

131.7218 

MOLS/KG. 


ENTHALPY OF THE 

MIXTURE 

• 144573E+09 

JOULES/KG. 


ENTROPY OF THE MIXTURE = 

. 332118E+Q5 

JOULES/KG. /KELVIN 

INT. ENERGY OF THE MIXTURE = 

» 1 33621 E+Q 9 

JOULES/KG. 


PRESSURE 

= 

. 109513E-03 

(BAR) 


RESI DUES : - 





0 . 

. 71054E-14 

0 . 0 

a 

0 . 

. 56843E-1 3 

. 22737E-12 

. 58797E-12 0 

• 

- . 35527E-13 


-09 

. 20269E-G2 
NO+ 

. 16440E-05 
. 13623E+0 7 
. 1 2792E+Q7 
.32018E+G3 
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Table II. Coefficients of Approximating Polynomials for Thermodynamic Functions 

[T ref = 0 K] 


TEMP. 

PANGE 

*(i,l) 

a : i,2 : 

a: ;,3 ■ 

a ; i ,4 ) 

• ( i ,5 ) 

•(1,6) 

a<i,7> 




1*1; Neutra i 

Diatonic Oxygen (02) 



AH 0 ^ 1 = 0 J/mol 


200. 

_ 

800. 

. 37703' 7 E*01 

- . 289522E-02 

. 953322E-05 

-.924699E-08 

.301919E-11 

- . 1885986+92 

.3693356+01 

800. 

- 

3000. 

. 289692E+01 

.237365E-02 

- . 149171E-05 

.4E6034E-09 

-.5394526-13 

. 8 22 4 6 4E +02 

.7501946+01 

3000. 

- 

6000. 

. 28*21 1E+01 

. 133206E-02 

- . 339159E-06 

.446522E-10 

-.2291486-14 

•5B3504E+03 

.862S50C+01 

6000. 

- 

10000. 

.S682UE+01 

- . 753006E-03 

. 229801E-06 

- .239559E-10 

.8004856-15 

-.2470046+44 

- .9874016+01 

10000. 

- 

15000. 

- .272590E+00 

.201151E-02 

-.245477E-06 

. 120252E-10 

-.2138996-15 

. 7611 77E +04 

. 3163176 *02 




i*2j Neutral 

Diatonic Nitrogen (N2> 



AH° f2 =0 J/mol 


200. 

_ 

806. 

. 346226E+01 

.582024E-03 

- . 305255E-05 

.622801E-08 

-.3375606-11 

. 8795176+00 

.3219266+01 

800. 

- 

3000. 

.270224E+01 

. 194439E-02 

- . 893000E-06 

. 197392E-09 

- . 1696786-13 

.2018026+03 

.7204086+01 

3000. 

- 

6000. 

.391435E+01 

.315371E-03 

- . 564810E-07 

. 360 123E-1 1 

.5235946-16 

- .5355206+03 

.2167206-01 

6000. 

- 

10000. 

. 126575E+01 

. 182698E-62 

- . 375839E-06 

.330337E-10 

- . 93651 16-15 

.3142706+04 

.1794336+02 

10000. 

- 

15000. 

.277159E+02 

-.7417356-62 

. 823960E-06 

* . 352853E-10 

. 4967126-15 

-.5694276+05 

-.1740296+03 




i*3; Neutral 

flonatomc Oxygen (0) 



AH° f3 =246 783 J/mol 

£00 . 

- 

800. 

.321671E+01 

- .378227E-02 

.847468E-05 

-.886582E-08 

. 353656E-11 

. 2964056+05 

.1852646+01 

S00. 

- 

3000. 

.260454E+01 

- . 172355E-03 

. 1 15741E-06 

-.364179E-10 

.460115E-14 

. 2972906+05 

.4586536+01 

3000. 

- 

6000. 

.281017E+01 

-.2903996-03 

. 908333E-07 

* . 994278E-1 1 

. 3770446-15 

. 29S3666+05 

. 32S365E+01 

6000. 

- 

10000. 

. 192093E+01 

. 21 7766E-03 

-.1828846-07 

.490510E-12 

.2850736-17 

. 3078346+05 

♦927486E+01 

10000. 

- 

15000. 

. 192093E+-01 

, 2177666-03 

-.1828846-07 

.4905106-12 

.2850736-17 

. 307834E+05 

•927486E+01 




i*4; Neutral 

Nitric Oxide 

(NO) 



AH 0 ,. =90 671 J/mol 
f,4 

200. 

_ 

800. 

. 420644E+01 

-.450984E-02 

. 105574E-04 

-.859194E-08 

.240471E-11 

♦ 1088906+05 

.2313796+01 

800. 

- 

3000. 

.275438E+01 

♦2309336-02 

- . 128234E-05 

.340435E-09 

-.3480756-13 

. 111343E+05 

.907897E+01 

3000. 

- 

6000. 

. 380 154E+01 

. 4985756-03 

125313E-06 

.140939E-10 

- . 4482026-15 

. 106666E+05 

.316194E+01 

6600. 

- 

10000. 

. 491 332E+01 

.617553E-06 

- . 552224E-07 

. 1 16865E-10 

-.546425E-15 

. 8845386+04 

-.4578696+01 

10000. 

- 

15000. 

.204567E+02 

- .614981E-02 

. 869143E-06 

-.508775E-10 

.106242E-14 

-.2295526+05 

-.1156206+03 




i *5; Neutral 

Monatoaic Nitrogen (N) 



AH° f5 =470 81 8 J/mol 

260. 


866. 

.2500006+01 

0. 

0. 

e. 

0. 

. 5662676+05 

.4180736+01 

800. 

- 

*>AAA 

JVW • 

.250751E+01 

-.2479796-04 

•29641SE-07 

-.1528816-10 

.2891376-14 

.5662496+95 

.4143196+01 

3000. 

- 

CAM 
DWV • 

.26376 IE +01 

-.8737336-05 

-.647727E-07 

.2347346-19 

-.1739626-14 

.5645236+05 

.3223216+01 

6000. 

- 

10000. 

.3372066+01 

-.885546E-03 

. 2529336-06 

- .231879E-10 

.7047146-15 

. 5627026+05 

-.1056406+01 

10000. 


15000. 

-. 1020566+02 

• 4293136-02 

-.493107E-06 

.2495116-10 

-.469699E-15 

. 849338E+05 

.9640376+02 




1*6; UnlpositWe Nitric Oxide (10+) 



AH° f6 =984 611 J/mol 

200. 


800. 

. 2565936+31 

- . 2' 7 4203E-03 

-.5228436-06 

. 302282E-08 

-.191127E-U 

. 118414E+06 

.3S93S7E+01 

800. 

- 

3000. 

.c?1526E+3l 

. 1899816-02 

- . 8534626-06 

. 1 839036-09 

- • 1532856-13 

. 118621E+06 

.7936916+01 

3000. 

- 

6000. 

♦ 3931106+01 

. 306608E-03 

-.5670S9E-07 

.455735E-11 

-.7776206-16 

.11786% +00 

•701982E+00 

6000. 

- 

10000. 

. 1065626+01 

.179180E-02 

-.330778E-06 

. 2505276-10 

-.S4S401E-1S 

.1221426+06 

.2032836+02 

10000. 

- 

15000. 

. 313^7 06 + 02 

- . 78l70" , E-02 

.754371E-06 

-.2473416-10 

.1595246-15 

.4812456+05 

-.20218S6+03 
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Table II. Concluded 


TEHP. RANGE •<».!) *(1,3) .<l.3> 


*•7; Electron (E-3 


200. 

- 800. 

.2500006401 

0. 

0. 

890. 

- 3000. 

.2590606+01 

0. 

0. 

3000. 

- 6000. 

.2500606401 

0. 

0. 

6000. 

- 10000. 

.2500006^01 

0. 

0. 

10000. 

- 15000. 

.2500006^01 

0. 

0. 



■Ci,4) 

• ' i,5 > e<i,6) 

• ( i,“M 



AH = 0 J/mol 


0. 

0. 

-.4249006-01 

-. 1173396402 

0. 

0. 

-.4249906-01 

-.1173396402 

0. 

0. 

-.4249906-01 

-.1173396402 

0. 

0. 

-.4249906-01 

-.1173396402 

0. 

0. 

-.4249906-01 

-.1173396402 





1*8; Jnipoaitlve flonatoeic Nitrogen (N+) 

A 

H° f8 =1 873 124 J/mol 

200. 

_ 

800. 

.3081706401 

-.3773756-02 

.9677156-05 

1 10995E-07 

•472486E-1 1 

.2253286406 

.2431596401 

800. 

- 

3000. 

.2534706401 

-.5506746-04 

.3644006-07 

-. 123132E-10 

. 192247E-14 

.2253956406 

.4789566401 

3600. 

- 

6000. 

.280056E401 

-.2956646-03 

.9664976-07 

-.1105366-10 

.4414526-15 

.2251856406 

.3130536401 

6000. 

- 

16000. 

.200477E401 

.1918976-03 

-.1549916-07 

.4116596-12 

.2628766-17 

.2262266406 

.8463506401 

10000. 

- 

15000. 

.2059846401 

. 1887546-03 

-.1807S5E-07 

•802843E-12 

-.1308966-16 

.2260286406 

.8025446401 




1-9; Uni posit wo flonatoeic Oxygen (04) 


AH° fg = 1 560 726 J/mol 

200. 


800. 

.2500006401 

0. 

0. 

0. 

0. 

.1877136406 

.4380206401 

800. 

- 

3000. 

.2501346401 

-.4043136-05 

.4327156-08 

-. 195277E-1 1 

. 316328E-15 

.1877136406 

.4373426401 

3000. 

- 

6600. 

.2225896401 

. 296501 E-03 

-.115611E-06 

.1846126-10 

-.892492E-15 

.1879126406 

.60S761E401 

6000. 

- 

10000. 

.4121106401 

-.867862E-03 

.1533756-06 

-.922601E-11 

. 178420E-15 

. 1854376406 

-.6638896401 

10000. 

- 

15600. 

.1734956401 

-.807107E-04 

.5980636-07 

-.4S9985E-I1 

. 102952E-15 

.1910746406 

.1079196402 




1-10; Neutral Argon (Ar) 



A H° f 1 Q =0 J/mol 


200. 

_ 

800. 

.2500606401 

6. 

0. 

0. 

0. 

-.4249906-01 

.4366506401 

800. 

- 

3000. 

.2560006401 

0. 

0. 

0. 

0. 

-.4249906-01 

.4366506401 

3000. 

- 

6000. 

.2500006401 

0. 

0. 

0. 

0. 

-.4249906-01 

.4366506401 

6000. 

- 

10000. 

.2630356401 

-.773934E-04 

.1721976-07 

170282E-11 

.6319886-16 

-.1754806403 

.3489016401 

10000. 

- 

15000. 

-.914240E401 

. 385082 E -02 

-.4680046-06 

.2450696-10 

-.4564696-15 

.2775196405 

.8946186402 




1*11; Uni posit Wo Argon (Ar4) 


AH° { 11 = 1 520 581 J/mol 

200. 

- 

800. 

.2681886401 

-.2072036-02 

.7478606-05 

-.8767206-08 

.3448426-11 

.182871E406 

.5076466401 

800. 

- 

3600. 

.2422066401 

.8634706-03 

-.7845946-06 

.2667466-09 

-.3185066-13 

.1828536406 

.S92363E401 

3000. 

- 

6000. 

.2558846401 

.8973626-04 

-.5009166-07 

.8829826-11 

-.5237196-15 

.1830146406 

.5531306401 

6000. 

- 

10000. 

.2501426401 

.2944666-04 

-.8047176-08 

.775358E-12 

-.2578096-16 

.1832526406 

.6054406401 

10000. 

- 

15000. 

.2570616401 

-.2062116-04 

.2978606-08 

-.2110196-12 

.5747686-17 

.1832236406 

.5616506401 
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Figure 1. Regions of applicability of various models. 




Figure 2. Log o^jo as function of temperature. 
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Figure 2. Concluded. 
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